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STRUCTURAL INSTABILITY OF NONLINEAR PLATES 
MODELLING SUSPENSION BRIDGES: 
MATHEMATICAL ANSWERS TO SOME LONG-STANDING QUESTIONS 

ELVISE BERCHIO, ALBERTO FERRERO, AND FILIPPO GAZZOLA 


Abstract. We model the roadway of a suspension bridge as a thin rectangular plate and we study in 
detail its oscillating modes. The plate is assumed to be hinged on its short edges and free on its long 
edges. Two different kinds of oscillating modes are found: longitudinal modes and torsional modes. 
Then we analyze a fourth order hyperbolic equation describing the dynamics of the bridge. In order 
to emphasize the structural behavior we consider an isolated equation with no forcing and damping. 
Due to the nonlinear behavior of the cables and hangers, a structural instability appears. With a finite 
dimensional approximation we prove that the system remains stable at low energies while numerical 
results show that for larger energies the system becomes unstable. We analyze the energy thresholds 
of instability and we show that the model allows to give answers to several questions left open by the 
Tacoma collapse in 1940. 


1. Introduction 

The history of suspension bridges essentially starts a couple of centuries ago. The first modern 
suspension bridge is considered to be the Jacob Creek Bridge, built in Pennsylvania in 1801 and 
designed by the Irish judge and engineer James Finley, see m for the patent and the original design. 
At the same time, several suspension bridges were erected in the UK, see (e.g.) the introduction in the 
seminal book m- The political instability due to the French Revolution and to the Napoleon period 
kept France slightly delayed. For this reason, M. Becquey (Conseiller d’Etat, Directeur General des 
Ponts et Chaussees et des Mines) committed Navier to visit the main bridges in the UK and to report 
on their feasibility and performances. In his detailed report [251 p.161], Navier wondered about the 
possible negative effects of the action of the wind: Les accidens qui resulteraient de cette action ne 
peuvent etre apprecies et prevenus que d’apres les lumieres fournies par Vobservation et Vexperience. 
Unfortunately, he had seen right. 

Many bridges manifested aerodynamic instability and uncontrolled oscillations leading to collapses, 
see e.g. hi nn. These accidents are due to many different causes and in this paper we are only 
interested about those due to wide unexpected oscillations. We will give a mathematical explanation 
for the appearance of torsional oscillations by analyzing a suitable partial differential equation modeling 
the bridge. 

Thanks to the videos available on the web |33j most people have seen the spectacular collapse of 
the Tacoma Narrows Bridge (TNB), occurred in 1940. In Figure [I] (picture taken from [2j P-6]) one 
can see the roadway of the TNB under a torsional oscillation. This kind of oscillation was considered 
the main cause of the collapse mm- But the appearance of torsional oscillations is not an isolated 
event occurred only at the TNB. The Brighton Chain Pier was erected in 1823 and collapsed in 1836: 
Reid [2B] reported valuable observations and sketched a picture illustrating the collapse see Figure [2] 
(picture taken from |25j). The Wheeling Suspension Bridge was erected in West Virginia in 1849 and 
collapsed in a violent windstorm in 1854: according to |34| . it twisted and writhed, and was dashed 
almost bottom upward. At last there seemed to be a determined twist along the entire span, about one 
half of the flooring being nearly reversed, and down went the immense structure from its dizzy height 
to the stream below, with an appalling crash and roar. Finally, let us mention that Irvine m Example 
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Figure 1. The collapsed Tacoma Narrows Bridge (1940). 



Figure 2. Collapse of the Brighton Chain Pier (1836). 


4.6, p.180] describes the collapse of the Matukituki Suspension Footbridge in New Zealand (occurred 
in 1977, just twelve days after completion) by writing that the deck persisted in lurching and twisting 
wildly until failure occurred, and for part of the time a node was noticeable at midspan. This description 
is completely similar to what can be seen in Figures [l] and [2] as well as in the video |33]. These are just 
few examples aiming to show that the very same instability was observed in several different bridges. 

These accidents raised some fundamental questions of deep interest also for mathematicians. Longi¬ 
tudinal oscillations are to be expected in suspension bridges but 

(Ql) why do longitudinal oscillations suddenly transform into torsional oscillations? 
This question has drawn the attention of both mathematicians and engineers but, so far, no unanimously 
accepted response has been found. The distinguished civil and aeronautical engineer Robert Scanlan 
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[29, p.209] attributes the appearance of torsional oscillations to some fortuitous condition. The word 
“fortuitous” highlights a lack of rigorous explanations and, according to m, no real progress has been 
done in subsequent years. 

The above collapses also show that the torsional oscillation has a particular shape, with a node at 
midspan. And it seems that this particular kind of torsional oscillation is the only one ever seen in 
suspension bridges. From the Official Report [2j p.31] we quote Prior to 10:00 A.M. on the day of the 
failure, there were no recorded instances of the oscillations being otherwise than the two cables in phase 
and with no torsional motions whereas from Smith-Vincent |32l p.21] we quote the only torsional mode 
which developed under wind action on the bridge or on the model is that with a single node at the center 
of the main span. This raises a further natural question: 

(Q2) why do torsional oscillations appear with a node at midspan? 

According to Eldridge [2] V-3], a witness on the day of the TNB collapse, the bridge appeared to 
be behaving in the customary manner and the motions were considerably less than had occurred many 
times before. From [2] p.20] we also learn that in the months prior to the collapse one principal mode 
of oscillation prevailed and that the modes of oscillation frequently changed. In particular, Farquharson 
[21 V-10] witnessed the collapse and wrote that the motions, which a moment before had involved a 
number of waves (nine or ten) had shifted almost instantly to two. This raises a third natural question: 

(Q3) are there longitudinal oscillations which are more prone to generate torsional 

oscillations? 

The purpose of this paper is to use the semilinear plate model developed in ESI and to adapt it to a 
suspension bridge having the same parameters as the collapsed TNB. By analyzing both theoretically 
and numerically this model, we will give an answer to the above questions (Ql), (Q2) and (Q3). 

This paper is organized as follows. In Section [2] we recall and slightly modify the model introduced 
in PCS]; in particular we discuss the nonlinear restoring force due to the hangers+cables system. In 
Section [3] we study in great detail the oscillating modes of the plate, according to the TNB parameters. 
In Section [4] we analyze the full evolution equation in the case where the system is isolated: we obtain 
a fourth order hyperbolic equations and we show that the initial-boundary-value problem is well posed. 
Then we define what we mean by torsional stability and we state two sufficient conditions for the 
stability. In Section [5] we numerically compute the thresholds of stability, according to our definition. 
In Section [6] we validate our results from several points of view: we show that the linearization and 
the uncoupling procedures do not alter the results and that a full numerical analysis does not give 
significantly different responses. Sections [7] and [8] are devoted to the proofs of the stability results. 
Finally, in Section [9] we afford an answer to the above questions. 

2. A NONLINEAR MODEL FOR A DYNAMIC SUSPENSION BRIDGE 

We follow the mathematical model suggested in [T6] by modifying it in some aspects. We view the 
roadway (or deck) of a suspension bridge as a long narrow rectangular thin plate hinged at the two 
opposite short edges and free on the remaining two edges. Let L denote its length and 2 1 denote its 
width; a realistic assumption is that 2£ = The rectangular plate H C l 2 is then 

n = (0, L) x (-£,£). 

Let us first discuss different positions of the plate depending on the forces acting on it. If the plate had 
no mass (as a sheet of paper) and there were no loads acting on the plate, it would take the horizontal 
equilibrium position uq, see Figure |3j If the plate was only subject to its own weight w (dead load) 
it would take a U-position such as u w in Figure [3] If the plate had no weight but it was subject to 
the restoring force of the cables-hangers system, it would take a lb-position such as u/,: this is also the 
position of the lower endpoints of the hangers before the roadway is installed. If both the weight and 
the action of the hangers are considered, the two effects cancel and the equilibrium position uq = 0 is 
recovered. 

Since the bending energy of the plate vanishes when it is in position uq = 0, the unknown function 
should be the displacement of the plate with respect to the equilibrium no- Augusti-Sepe [5] (see also 
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Figure 3. Different positions of the bridge. 


0 ) view the restoring force at the endpoints of a cross-section of the roadway as composed by two 
connected springs, the top one representing the action of the sustaining cable and the bottom one 
(connected with the roadway) representing the hangers, see Figure |4j 


< cables-> 

< -hangers-> 


Figure 4. The cables+hangers system modeled with two connected springs. 

The action of the cables is considered by Bartoli-Spinelli j6j p. 180] the main cause of the nonlinearity 
of the restoring force: they suggest quadratic and cubic perturbations of a linear behavior. If u denotes 
the downwards displacement of the deck, here we simply take 

(1) g(u) = ki u + &2 u 3 

for some k\,k 2 >0 depending on the elasticity of the cables and hangers. Let us mention that Plaut- 
Davis pT i, § 3.5] make the same choice. 

The action of the hangers on the roadway is confined in the union of two thin strips parallel and 
adjacent to the two long edges of the plate D, that is, in a set of the type 

(2) lj:={0,L)x[(-Z,-Z + s)U(Z-£,Z)] 

with e > 0 small compared to Z. Summarizing, we take as restoring force and potential due to the 
cables-hangers system 

(3) h(y,u) = T(y) {k 1 u + k 2 u^j , H(y,u) = h(y, r)dr = T(y) u 2 + y , 

where T is the characteristic function of {—Z, — Z + e) U (Z — £,Z). 

The derivation of the bending energy of an elastic plate goes back to Kirchhoff |22] and Love [2T] , 
see also [16] for a synthesized form. The total static energy of the bridge is obtained by adding the 
potential energy from ([3]) to the bending energy of the plate: 

(4) E t (u) = u ^ d _ ^ J ^ Q (Au) 2 + (<r - 1) det(ZAi)^ + (H(y,u ) - fu ) 

where / is an external force, d denotes the thickness of the plate, E is the Young modulus and a is the 
Poisson ratio. For a plate, the Poisson ratio is the negative ratio of transverse to axial strain: when a 
material is compressed in one direction, it tends to expand in the other two directions. The Poisson 
ratio a is a measure of this effect, it is the fraction of expansion divided by the fraction of compression 
for small values of these changes. Usually one has 

o<„<L 


( 5 ) 
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A functional space where the energy E t is well-defined is 

H 2 (Q) : = {w € w = (1 on {0 , L} X (-£,£)} 

which is a Hilbert space when endowed with the scalar product 

(6) (it, v ) H 2 := / [AuAv + (1 - a)(2u xy v xy - u xx v yy - u yy v xx )] dxdy , 

J n 

see m Lemma 4.1]. We also consider 

%{Sl) := the dual space of H 2 (fl) 

and we denote by (•, •) the corresponding duality. Since we are in the plane, H 2 (fl) C C°(fi) so that the 
condition on {0,L} x (—£,£) introduced in the definition of H 2 (Q) is satisfied pointwise. If / G L 1 (H) 
then the functional E t is well-defined in H 2 (Q), while if / G %{£t) we need to replace fu with (/, u) 
although we will not mention this in the sequel. 

If the load also depends on time, / = f(x,y,t), and if m denotes the mass density of the plate, then 
the kinetic energy of the plate should be added to the static energy 

(7) -.. ' * E6? 


„ . , m f 

Su(t) ■■= y 


u\ + 


+ (o- - 1) det(L> 2 «)^ + ( H (y, u ) - fu) 


/n 12(1 — a 2 ) 

This is the total energy of a nonlinear dynamic bridge. As for the action, one has to take the difference 
between kinetic energy and potential energy and integrate over an interval of time [0, T\: 


rT r 


v4.(it) : = 


Ed 3 


m / 2 

2 J n Ut 12(1 — a 2 ) 


(A u) 


+ (cr - 1) det {D 2 u) ) - (; y , u ) - fu) 


dt. 


The equation of the motion of the bridge is obtained by taking the critical points of the functional A: 

Ed 3 


m u t t + 


A~u + h(y, u) = f in Q x (0, T). 


12(1 -a 2 ) 

Due to internal friction, we add a damping term and obtain 

mutt+5u t + l2 ^ a2) A 2 u+h(y,u) = f in fix(0,T) 

it(0, y, t) = u xx { 0, y, t) = u(L , y, t)=u xx (L, y,t) = 0 for (y, t) G (-£, £) x (0, T) 
u yy (x, ±£, t)+au xx (x, ±£, t) = 0 for (.x , t) G (0, L) x (0, T) 


( 8 ) 


u. 


yyy 


(x, ±£,t) + ( 2 - a)u xxy (x,±£,t) = 0 


u(x,y,0)=u 0 (x,y ), u t (x,y,0) = u 1 (x,y) 


for (x, t) G (0, L) x (0, T) 
for (x, y) G D 


where 5 is a positive constant. We refer to HS1 for the derivation of the boundary conditions. 

For a different model of suspension bridges, similar to the one considered in [3], Irvine [19] p.176] 
ignores damping of both structural and aerodynamic origin. His purpose is to simplify as much as 
possible the model by maintaining its essence, that is, the conceptual design of bridges. Here we follow 
this suggestion and consider the isolated version of ([8]) for which global existence is expected (T = oo). 
This isolated version of (l8|) reads 


Ed 3 


-A 2 w+T(y)(kiw + k 2 iv 3 ) = 0 for (x, y) G (0, L) x (—£, £) , t > 0 . 


(9) ™“«+ 12(w2) 

where w denotes the downwards vertical displacement and all the constants are defined in Section [2j 
In order to set up a reliable model, we consider the lengths of the plate as in the collapsed TNB. 
According to |2j p.ll], we have 

21 _ 39 ^ 1 _ ^ 

L 2800 ~ 75 vr ' 


(10) L = 2800 ft. « 853.44 m, 2£ = 39 ft. » 11.89 m , that is, 
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Therefore, we may scale the plate (0,L) x (—£,£) to (0,7r) x (—jfg, -jf^). Then we take the amplitude 
of the strip oj (see ([ 2 ])) containing the hangers also as at the TNB, see [2j p. 11]: 

. . 7T 

( 11 ) £= 


1500 


Referring to Section [5] for the values of the involved parameters, we put 


( 12 ) 


iki / ttx ny ki 
w[x, y, t) = \ — u\ 


L L 


and 


Ed 3 


7 r 


m 


7 = 


\2k\(l — a 2 ) L 4 ' 


Then ([ 9 ]) becomes an equation where the only parameter is the coefficient of the biharmonic term: 


u tt + 7 A 2 u + T(y) (u + u 3 ) = 0 in (0, ir) x 


7T 7T \ 

150’150/ 


x 


and T is the characteristic function of the set (— ifo>~Mi) U (iif)’ifo) - Finally, we notice that for 
metals the value of a lies around 0.3, see [25 p. 105], while for concrete we have 0.1 < a < 0.2. Since 
the suspended structure of the TNB consisted of a “mixture” of concrete and metal (see [2 P- 13]), we 
take 


(13) 


a = 0 . 2 . 


(14) 


Then by using (12) we find the dimensionless version of the problem under study 

u tt +^A 2 u+T(y)(u + u 3 ) = 0 in 17 x (0, 00 ) 

u(0,y,t) = u xx (0,y,t)=u(Tr,y,t)=u xx (ir,y,t)=0 for (y,t)e (-^ I f 5 )x(0,oo) 
u yy (x, ± I | g ,f)+0.2 • u xx (x,±j^,t) = 0 for (x, t) £ (0, n) x (0, 00 ) 


l yyyy 


(x, ± 15Q , t) + 1.8 • U X xy (37 A f ) 0 


u(x,y,0)=u o {x,y ), u t (x,y,0) = u 1 (x,y) 


for (x, t) G (0,7r) x (0, 00 ) 
for (x, y) G 


where = (0,7r) x (—xfo’ lfo)' The initial-boundary value problem (|14h is isolated, which means that 
it has a conserved quantity. This quantity is the energy introduced irijjTj) which is constant in time: 

(15) E(u) = ^u\ dxdy + Q(Au) 2 + y {u 2 xy - u xx u y y ) + T(y) ("T + t)) dxdy ' 

3. The eigenfunctions of the linearized problem 

For the rectangular plate = (0, n) x (—£,£) with Poisson ratio a we are interested in the eigenfunc¬ 
tions of the eigenvalue problem 


(16) 


A 2 w = X w 

w(0, y) = Wxxitt, y) = w(tt , y) = w xx { tt, y) = 0 

Wyy(x, ±1) + aw X x(x,±£) = w yyy (x ,±£) + (2 - a)w XX y(x, ±£) = 0 

Problem (16) admits the following variational formulation: a nontrivial function w E is an 

eigenfunction of (|16[) if there exists A E M (an eigenvalue) such that 


in 12 

for y E {—£,£) 
for x E (0,7r). 


/ [AwAv + (1 — a){2w xy v X y — w xx v y y — Wy y v xx ) — A wv] dxdy = 0 for all v E iP 2 (fl) . 

•/>> 

We recall from m Theorem 7.6] a statement describing the whole spectrum and characterizing the 
eigenfunctions. It is shown there that the eigenfunctions may have one of the following forms: 

Proposition 3.1. Assume <©■ Then the set of eigenvalues of ( fib] ) may be ordered in an increasing 
sequence {A*,} of strictly positive numbers diverging to +oo and any eigenfunction belongs to C°°(Q); 
the set of eigenfunctions of (16) is a complete system in Moreover: 
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(i) for any m ^ 1, there exists a unique eigenvalue A = i G ((1 — cr) 2 m 4 , ?n 4 ) with corresponding 
eigenfunction 


r 1/2 n , 2i cosh {W m 2 + ^ i) , r 1/2 , „ ^ 2i cosh [y^-ulZ 

Ki - (! - J - ) 7 „ 1/0 < + Kl,i + (! - J 


cosh 




cosh I ^ 




sin(m.T); 


(ii) /or any m ^ 1, f/iere exist infinitely many eigenvalues A = ti m ,k > m 4 (k ^ 2) with corresponding 
eigenfunctions 


[fill - (! - (j ) m2 ] 


cosh [y\]pj^k+rrt 


+ + (! “ CJ ) m2 ] 


cosh I € 


VS 


+m z 


COS I 

[y\/»™ 2 k -™ 2 ) 

cos 1 

(V^fc-™ 2 ) 


sin(mx); 


(Hi) for any rri ^ 1, 1/iere exist infinitely many eigenvalues A = u m _k > m 4 (k ^ 2) with corresponding 
eigenfunctions 


["til ~ C 1 - <r)™ 2 ] 


sinh (yy/ v ^k +m2 ) , r ,1/2 ^ n _ ^_21 sin (W 


sinh t 


i\[7^ 

V m > 


+m z 


+ u - ct )™ 2 ] 


1/2 o 
z/' u—'m 


(V- 


iW“ — m 2 

m.K 


sm(mx) 


(iv) for any m ^ 1 satisfying tm\J2 coth(£m\/ 2 ) > (A ^) 2 tZiere exists an eigenvalue A = z>Vn,i G 
, m 4 ) with corresponding eigenfunction 


r 1/2 n ^ 21 sinh (V m 2 +^i) , r 1/2 , n , 21 Sinh (W m2 -' / ™.i 

L<i - (! - CT )™ J - ) , 1/9 \ + K ,1 + (! - CT )™ J 


sinh I f-i/ m 2 +ulf 2 . 


sinh I f.\!m 2 —vlf 2 . 


sin(mx). 


Finally, if the unique positive solution s > 0 of the equation 


(17) 


tanh(\/ 2 s£) = 


(7 


2-o- 


\/2s£ 


is not an integer, then the only eigenvalues and eigenfunctions are the ones given in (i) — (iv). 

Of course, © has probability 0 to occur in a real bridge; if it occurs, there is an additional eigenvalue 
and eigenfunction, see |16j . The eigenvalues Afc are solutions of explicit equations. More precisely: 

(i) the eigenvalue A = is the unique value A G ((1 — <x) 2 m 4 , m 4 ) such that 


\J m 2 —A 1 / 2 (A 1 / 2 + (l —o-)m 2 ) 2 tanh(£\/ m 2 — A 1 / 2 ) = \/ m 2 + A 1 / 2 (A 1 / 2 — (1 —er)m 2 ) 2 tanh(£\/ m 2 + AV2); 
(ii) the eigenvalues A = p m ,k (k ^ 2) are the solutions A > m 4 of the equation 
\J A 1 / 2 — m 2 (A 1 ^ 2 + (l— a)m 2 ) 2 t&n(£\/ A 1 / 2 — m 2 ) = —\J A 1 / 2 +m 2 ( A l/2_ (1 _ a) 

m 2 ) 2 tanh(£V / AV2 + m 2 ); 

(Hi) the eigenvalues A = v m ,k (k ^ 2) are the solutions A > m 4 of the equation 
\J A 1 / 2 — m 2 (A 1//2 + (l — a)m 2 y tanh(£-\/A 1 / 2 +m 2 ) = VA 1 / 2 +m 2 (A 1 / 2 — (1 — a)m 2 ) 2 tan(^-\/A 1 / 2 —m 2 ); 
(iu) the eigenvalue A = z/ m .i is the unique value A G ((1 — cr) 2 m , m 4 ) such that 

\Jm 2 - A 1 / 2 (A 1 / 2 + (1 - o-)m 2 ) 2 tanh(^\/ A 1 / 2 +m 2 ) = \/ A 1 / 2 +m 2 (A 1 / 2 — (1 —er)m 2 ) 2 tanh(^V / A 1 / 2 — m 2 ) . 

The least eigenvalue is Ai = pi.i- the corresponding eigenfunction is of one sign over and this fact 
is by far nontrivial. It is well-known that the first eigenfunction of some biharmonic problems may 
change sign. When 11 is a square, Coffman m proved that the first eigenfunction of the clamped plate 
problem changes sign, see also [[21] for more general results. Moreover, Knightly-Sather 123 Section 
3] show that the buckling eigenvalue problem for a fully hinged (simply supported) rectangular plate, 
that is with u = A u = 0 on the four edges, may admit a least eigenvalue of multiplicity 2. Hence, the 
positivity of the first eigenfunction of (16) is not for free. Due to the L 2 -orthogonality of eigenfunctions, 


it is the only positive eigenfunction of (16) 
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The eigenfunctions in (i) — (ii) are even with respect to y whereas the eigenfunctions in (Hi) — (iv) 
are odd. We call longitudinal eigenfunctions the eigenfunctions of the kind (i) — ( ii ) and torsional 
eigenfunctions the eigenfunctions of the kind (in) — (iv). Since i is small, the former are essentially 
of the kind c rn sin(rrtx) whereas the latter are of the kind c m ysm(mx). The pictures in Figure [5] 
display the first two longitudinal eigenfunctions (approximately described by c\ sin(x) and C2sin(2x)) 
and the second torsional eigenfunction (approximately described by C2ysin(2x)). In each picture the 
displacement of the roadway is compared with equilibrium. 



Figure 5. Some possible oscillations of a bridge roadway. 


Of particular interest is the lower picture in Figure [5] which corresponds to a torsional eigenfunction 
with a node at midspan. This is precisely the behavior of the TNB prior to its collapse in November 
1940, see the video |33] and Figure [lj This is also the behavior of the Brighton Chain Pier, see Figure 
[2j and of the other bridges described in the introduction. With the notations of Proposition 3.1, this 
common behavior of suspension bridges may be rephrased as follows 


the oscillations causing the collapse of a suspension bridge 
are of the kind C2?/sin(2x), as represented in the bottom picture of Figure [5j 
and correspond to the eigenvalue 1 / 2 , 2 - as given by Proposition |3.1| 


(18) 


We remark that the eigenfunction corresponding to the eigenvalue z/ 1,1 is of the kind C 2 ?/sin(x) but 
this eigenvalue in general does not exist since the inequality in Proposition 3.1 (iv) is usually satisfied 
only for large m. 


By (10) we may take 


(19) 


7r 

150 


With the choices in (13) and (19) we numerically obtained the eigenvalues of (16) as reported in Table 
[I] We only quote the least 16 eigenvalues because we are mainly interested in the second torsional 
eigenvalue which is, precisely, the 16th. 


Our results are obtained with the parameters of the TNB, see (10), (13) and (19). As already 


mentioned in the introduction, Farquharson j2j V-10] witnessed the collapse and wrote that the motions, 
which a moment before had involved a number of waves (nine or ten) had shifted almost instantly to 
two. Note that the longitudinal eigenvalue immediately preceding the least torsional eigenvalue is faio,i- 
it involves the function sin(10x) which has precisely “ten waves”. This explains why at the TNB the 
torsional instability occurred when the bridge was longitudinally oscillating like sin(10x). Then, if 
no constraint acts on the deck, the energy should transfer to the eigenfunction corresponding to 1 / 1,2 
which has a behavior like ysin(x). Among longitudinal eigenfunctions, the tenth is the most prone 
to torsional instability since the ratio between torsional eigenvalues and longitudinal eigenvalues is 
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eigenvalue 

Ai 

A2 

A3 

a 4 

A5 

Ae 

A7 

As 

kind 

AH,i 

A*2,l 

h3,l 

/M,l 

Mo,i 

W,i 


V 8,1 

Veigenvalue ~ 

0.98 

3.92 

8.82 

15.68 

24.5 

35.28 

48.02 

62.73 


eigenvalue 

Ag 

A 10 

An 

Al2 

A 13 

A 14 

A 15 

Al6 

kind 

fJ-9,1 

M 10,1 

2 

y-11,1 

M 12,1 

hl3,l 

Ml4,l 

^2,2 

V eigenvalue ~ 

79.39 

98.03 

104.61 

118.62 

141.19 

165.72 

192.21 

209.25 


Table 1. Approximate va' 


ue of t 


le least 16 eigenvalues of (16) for a = 0.2 and l 


7T 

150* 


minimal (close to 1) precisely for the tenth longitudinal eigenfunction: why small ratios yield strong 
instability is explained for a simplified model through Mathieu equations in [9]. However, in the case 
of a bridge, the sustaining cable yields a serious constraint. With a rude approximation, the cables 
may be considered as inextensible. A better point of view is that they are only “weakly extensible”, 
which means that their elongation cannot bee too large. In Figure [6] we represent the deformation 



Figure 6 . Elongation of the cable generated by the oscillations of the deck. 


of a cable in the two situations where the deck behaves like sin(x') and like sin(2x). It turns out 
that the no-noded behavior sin(x) (on the left) only allows small vertical displacements of the deck 
(between the grey and red positions) and, therefore, small torsional oscillations of the kind y sin(x). 
This is confirmed by numerical experiments for models where the cable plays a dominant role, see m 
where it is shown that, basically, “the first mode does not exist” in actual bridges. On the contrary, 
the one-noded behavior sin(2x) allows much larger torsional oscillations of the kind ysin(2x), see the 
right picture. Our explanation of the transition described by Farquharson is that when the oscillation 
sin(10x) became sufficiently large, reaching the threshold of the torsional instability, the cable forced 


the transition to the eigenfunction y sin(2x) instead of y sin(x). This gives a sound explanation to (18) 
and a first answer to (Q2), see Section [5] 

If we slightly modify the choices in (|13[ ) and (19), still in the range of the TNB, we obtain the 
eigenvalues of (16) as reported in Table |2| 


While comparing with Table [TJ it is noticeable that all the eigenvalues have slightly lowered but the 
qualitative behavior and the corresponding explanation remain the same. 


4. Torsional stability of the longitudinal modes 

4.1. Existence, uniqueness, and finite dimensional approximation of the solution. Let T be 

the characteristic function of the set (— 530 , — ^ 3 ) U ( 353 , 333 ) anc ^ ^ 


h(y,u) = T(y) (u + u 3 ) . 


( 20 ) 

We say that 
( 21 ) 


u G C°(M+;i/, 2 (L?))nC 1 (M + ;L 2 (II))nC 2 (M + ;H(II)) 













































10 


ELVISE BERCHIO, ALBERTO FERRERO, AND FILIPPO GAZZOLA 


eigenvalue 

Ai 

A 2 

A 3 

a 4 

A 5 

^6 

A 7 

^8 

kind 

Mi,i 

M2,l 

^3,1 

M4,l 

M5,l 

M6,l 

M 7,1 

M8,l 

Veigenvalue ~ 

0.97 

3.87 

8.71 

15.49 

24.21 

34.87 

47.46 

62 


eigenvalue 

Ag 

A 10 

An 

Al2 

A 13 

Al4 

A 15 

Al6 

kind 

^ 9,1 

Aho,i 

Pl ,2 

Mh,i 

^12,1 

Ml3,l 

Ml4,l 

P2,2 

Veigenvalue ~ 

78.48 

96.9 

97.24 

117.27 

139.58 

163.84 

190.1 

194.51 


Table 2. Approximate value of the least 16 eigenvalues of 


16) for a 


0.25 and l = 


7T 

144* 


is a solution of (14) if it satisfies the initial conditions and if 

(22) (u”(t),v) + 7 {u(t),v) H 2 + ( h(y , u(t)),v ) L 2 =0 \/v G H 2 (n ), Vf G (0, T ), 

see ©>• Then, by arguing as in m Theorem 3.6] we may prove the following result. 

Theorem 4.1. Assume ([5]). Let uq G and u\ G L 2 (Ll). Then there exists a unique solution 

u = u(t) of (14) and its energy (15) satisfies 

E (u(t)) = j ^ dxdy + ^(Au 0 ) 2 + y ((u 0 ) 2 xy - (u 0 ) xx (u 0 )yy) + T (y) (y + y)) dxd V ■ 

Sketch of the proof of Theorem f.l The proof of Theorem 4.1 makes use of a Galerkin method. The 
solution of (14) is the limit (in a suitable topology) of a sequence of solutions of approximated problems 
in finite dimensional spaces. By Proposition |3.1| we may consider an orthogonal complete system 


{w'fcjfc^i C of eigenfunctions of (16) such that ||uifc || L 2 = 1. Let {\k\k^i be the corresponding 

eigenvalues and, for any mfi 1, put W m := span{u;i,..., w m }. For any m 1 let 


u 0 : — 


:= ^( u 0 ,Wi) L 2Wi = 1 (u 0 ,w i ) H 2 Wi and ufi = Wi 


2=1 


2=1 


2=1 


so that u™ —> uq in and — > u\ in L 2 (fii) as m —>• + 00 . Fix T > 0; for any m 1 one seeks a 

solution u m G C 2 ([0,T]; W m ) of the variational problem 


(23) 


(u"(t),v) L 2 +'y(u(t),v) H 2 + (h(y,u{t)),v) L 2 = 0 
u{ 0 ) = < , u'( 0 ) = ufi 


for any v G W m and t G (0, T). If we put 

m 

(24) u m {t) = Y J 9 ?{t)w i and g m (t) := ..., g™(t)) T 


2=1 


then the vector valued function g m solves 

(25) I (gm(t)) " + 7 Amgm ' (t) + ® m(gmiW) = ° ^ G ( °’ T) 

j 5 m (0) = ((u 0 ,uii) L 2,...,(u 0 ,u; m ) L 2) T , (5 m )'(0) = ((ui,u;i) L 2,...,(ui,u; m ) L 2) T 

where A m := diag(Ai,..., A m ) and <l> m : M m — > M m is the map defined by 


I m m 

$rn(£l,---,£m) ■= 


3 = 1 


1=1 


L 2 
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From (20) we deduce that <F m e Lipi 0 


) and hence (25) admits a unique solution. Whence, the 


function u m (t) in (24) belongs C* 2 ([0, T); H 2 (Q)) and is a solution of the problem 


(26) 


+7 Lu m (t) + Pm(h(y,U m (t))) = 0 
u m ( 0 ) = ufi l , u' m ( 0 ) = u™ 


for any t ^ 0 


where L : H 2 (Q) — »• %{Sl) is implicitly defined by (Lu,v) := (u,v)h 2 for any u, v G H*(Q), and P m 
is the orthogonal projection from H 2 (fl) onto W m . By arguing as in [16], one finds that the sequence 

□ 


{um} converges in C'°([0, T\; H*(Q)) D C' 1 ([0, T\; L 2 (fi)) to a solution of (14). 

The above proof shows that the solution of (|14|) may be obtained as the limit of a finite dimensional 


analysis performed with a finite number of modes. Let us fix some value E > 0 for (15). Depending 


on the value of E, higher modes may be dropped, see [ 8 ], Section 3.3] for a detailed physical motivation 
and for some quantitative estimates. Our purpose is to study the torsional stability of the low modes 
in a sense that will be made precise in next section. In particular, the results obtained in Section [3] 
suggest to focus the attention on the lowest 16 modes, including the two least torsional modes. 

4.2. A theoretical characterization of torsional stability. In this section we give a precise defi¬ 
nition of torsional stability. We point out that this might not be the only possible definition. However, 
the numerical results reported in Sections [5] and [ 6 ] show that our characterization well describes the 
instability. 

We take again h (y, u ) as in (120) and we fix m = 16 which is the position of the second torsional mode. 
From Proposition 3.1 and Table ijwe know that the eigenvalues and the L 2 -normalized eigenfunctions 
up to the 16th are given by 


Afc — < 


fJ’k, 1 
PL,2 
Wfc-1,1 
^2,2 


if 1 < A: < 10 
if k = 11 
if 12 < k < 15 
if k = 16 


Wk(x,y) = < 


•Vk(y)sin(kx) 

81 (y) sinpr) 

H )1 

v k—l{y) sin((fc—l)x) 
u k -i 
d 2 (y) sin(2a:) 
u> 2 


if 1 ^ k ^ 10 
if k = 11 
if 12 O ^ 15 
if k = 16 


with 


Vkiy) ■= \-Pk 


cosh yjf3 


cosh 


+ 


« - TJ 


k 2 1 cosh [y<JP l 


cosh 


h{y) ■= 


k 2 


sinh r , k 2 ~\ sin [y\/ak ) 

-T- otl~— - 

^ sin 


sinh < iisv“fe 


(fc = 1,..., 14), 


(A: = 1,2). 


where j3 k := k 2 ± fi k l (for k = 1,..., 14), a.± := vj 2 ± k 2 (for k = 1,2) and 


(27) u% = ir v k (y) dy (k = 1,..., 14) , Jf. = tt / 0%(y) dy {k = 1,2). 

Jo Jo 

Notice that the v k are even with respect to y, while 9\ and 62 are odd. 


Following the Galerkin procedure described in the proof of Theorem 4.1, we seek solutions of (23) in 
the form 

. . v-' / , vu(y)sin(kx) . , 9i.(y)sm(kx) 

u(x,y,t) = J2Mt) Js ^ 1 -—+ K f 

• Mi. • 


fc=i 




k =1 




where the functions (p k and 77 are to be determined. Take h as in (20) and, for all ( 77 ,..., 9214 , n, 72 ) € 


>16 


, put 


$fc(<Pi, •••) 7*14) t~i,T 2 ) = (h(y, fj 


14 2 

%‘C») sin Ci# I ^ 8 j(y)sm(jx) ^ y k (y) sin(fcr) 


3 =1 


i=i 


3 LJj 


L 2 


{k = 1,..., 14) 
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14 


rfc(¥>i, -,^14 ,ti,t 2 ) = [h[y,^2,ipj 

i=i 


2 

i’j( y) sin(j.r) gjfa) sin(jx) 

OJj ' / V j UJj 

3 = 1 


^fc(y) sin(fcx) 




L 2 


(A; = 1,2). 


(28) 


Then (26) becomes the system of ODE’s: 

+7Mfe,i^fe(A) + $k{<Pi{t), ...,Lpu(t),Ti(t),T 2 (t)) =0 (k = 1,14) 

Tk(t) +1 Vk,2Tk{t) + r k (ip 1 (t),...,ipu(t),Ti(t),T 2 {t)) =0 (A; = 1,2) 

for all t E (0, T). For all 1 ^ k ^ 14 we put '&k( l Pk{t)) = 4>*,(0, <pk(t), ..., 0). By taking into account 
that 

7T f^ 37T 

/ sin 2 (/cx) dx = — , / sin 4 (/cx) dx = — , Vk ^ 1 , 

Jo 2 Jo 8 

and that Vk is even with respect to y, some computations yield 

(29) 
where 

(30) 


^k(<Pk(t)) = a k ipk(t) + b k (fl(t ), 


7r 


150 


/i _ t 37T 

ak = ^9 / v'k(y)dy and b k = — T 
ujr J An i 4a;? 

k J 500 « 


' 37T 
500 


vt(y) d y (k = i, •••> 14) 


In particular, by combining (30) with (27) we see that 


Vjfe 


&k — 


l 2 (^~ 

V500’150 > ^ i 




L 4(3w_ _s_\ 
1 500 ’ 150 > 


I^IH,2( 0>I |_) 47r ll^fc|l I ,2 (0i _|_ ) 

We may now define what we mean by longitudinal mode. We point out that this is a classical 
definition in a linear regime while it is by no means standard how to characterize modes in nonlinear 
regimes; contrary to the linear case, the frequency of a nonlinear mode depends on the energy or, 
equivalently, on the amplitude of its oscillations. 


Definition 4.2. Let 1 ^ k ^ 14, M 2 3 (</>O)0i) ^ (0,0) and 4 'k as in (29). We call k-th longitudinal 
mode at energy E((J)q, (f>\) > 0 the unique (periodic) solution Tp k of the Cauchy problem: 


(31) 


f ^l'(A) + + ^k(<Pk(t)) = 0 

\¥>fc(O) = 0§, <p k (0) = <t>\ ■ 


Vi > 0 


If it were = 0 then the equation (31) would be linear and Definition 4.2 would coincide with the 


usual one: in this case, there would be no need to emphasize the dependence on the energy since the 
solution with initial data (f k {0) = and ip' k ( 0 ) = a(j)\ (for any a) would coincide with the solution 
of (31) multiplied by a. In view of (29), we have instead a nonlinear equation and (31) becomes 


(32) 


f Vk{t) + (7MM + a-k)<Pk(t) + b k ipl(t) = 0 Vf > 0 
\<Pk( 0) = ^o, <Pfc(°) = ■ 


The system (32) admits the conserved quantity 
(33) 


E = Cf + (7WM +at fl +b jA^ = M ' 2 


(00 


k\ 2 


+ (7Mfc,i + Ofc)—~-k bk 


(00 


,fc\4 


Any couple of initial data having the same energy leads to the same solution of (32) up to a time 


translation while it is no longer true that multiplying the initial data by a constant leads to proportional 


solutions; it is well-known that different energies yield different frequencies of the solution, see also (53) 
below. 

In order to define the torsional stability of a longitudinal mode tp k , we linearize the last two equations 
of system (28) around (0, ...,ip k (t ),..., 0) E M 16 . These two equations correspond, respectively, to the 
first and second torsional mode. In both cases we obtain a Hill equation of the type 

(34) £"(*) + A l>k (t)£(t) = 0 , 
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where, for every 1 ^ k ^ 14 and Z = 1, 2, we set 

(35) A i}k (t) = 7^,2 + ai + di jk lpl(t) 


with 


7T / 150 0 

(36) a* = -2 / 0?{y)dy = 


r 2 / 37r 7r \ 
^ 500’150 

^ 2 (0.lfo) 


< 1 , = ‘ 


*57/if r$<,v)0?(»)dy if l = k 

l l 500 

*57 /»»!(»)«?(!/) <7 if ^ • 

fc ( 50Q 


For the computation of these coefficients we used the facts that the integrals containing odd powers of 
Oi{y) vanish and that 

2 2 i tF if l = k 

sut(Z.t) shr(Zcx) dx = 

' f m^k. 


/ 0 


Since (34) is a linear equation with periodic coefficients the notion of stabiiity of its trivial solution 


is standard. This enables us to define the torsional stability of a longitudinal mode. 

Definition 4.3. Fix 1 ^ k 14 and l = 1,2. We say that the k-th longitudinal mode Tp k at energy 
E(f> q, <(>\), namely the unique periodic solution of (32), is stable with respect to the Z-th torsional 
mode if the trivial solution of (34) is stable. 


4.3. Sufficient conditions for the torsional stability. It is well-known that the stability regions 
for the Hill equations may have strange shapes such as pockets and resonance tongues, see e.g. nans]- 
Therefore, the theoretical stability analysis of any such equation has to deal with these shapes and with 


the lack of a precise characterization of the stability regions. For (34), the theoretical obstruction is 
essentially related to the following condition 


(37) 


I 7 * 4,2 + ai 

7Wfc,i + a k 


N 


where 7 is defined in (12) while a k and ai are defined, respectively, in (30) and in (36). The usual 


difficulties are further increased for (34) which, instead of a single equation, represents a family of Hill 


equations having coefficients with periods depending on the energy of the solution of (32) 


Below we discuss in some detail assumption (37). But let us start the stability analysis with the 


following sufficient condition for the stability of a longitudinal mode. 


Theorem 4.4. Fix 1 ^ k ^ 14, Z £ {1,2} and assume that (37) holds. Then there exists E l k > 0 and a 
strictly increasing function A such that A(0) = 0 and such that the k-th longitudinal mode Tp k at energy 

E(f k 


0! 


t>i) (that is, the solution of ([32 ]) ) is stable with respect to the l-th torsional mode provided that 


E < E 


or, equivalently, provided that 


< A(4) 


Theorem 4.4 is not a perturbation result. The pr oof g iven in Section [7] allows us to determine explicit 

in a qualitative form in order not to spoil the 
for the precise value of E\. and A (E l k ). 


4.4 


values of El and A (El). Here we stated Theorem 
statement with too many constants. We refer to Theorem 


7.3 


Assumption (37) is a generic assumption, it has probability 1 to occur among all random choices of 
the positive real numbers 7 , 2 , ai, a k - But even in the case where (37) fails we may obtain a 

sufficient condition for the torsional stability of longitudinal modes. 

Theorem 4.5. Fix 1 ^ k ^ 14, Z £ {1,2} and assume that there exists m £ N such that 


(38) 


7*4,2 + ai 
7A*fc,i T a k 


= (m + l) 2 . 
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Assume moreover that 
(39) 


2^2 + (m + l)7r^ di t h < 37r(m + l) 3 6fc 


TTien f/ie same conclusions of Theorem 4-4 hold. 


Theorem 4.5 which is proved in Section [8] raises the attention on the further technical assumption 
(39). We are confident that it might be weakened (see Remark 8.4 at the end of Section[8]) and, perhaps, 
completely removed (see 0). The reason is that, although they are not explicitly known, the stability 
regions for the Hill equations are very precise and, for the model problem (|34|), there is an additional 


energy parameter which could vary the stability regions. However, we will not discuss (39) here. 


Theorems 4.4 and 4.5 state that a crucial role is played by the amount of energy inside the system. In 


the next sections we make some numerical experiments on the equations (34) and we study the stability 


of the least 14 longitudinal modes of the plate with the TNB parameters. Our results show that for 
each longitudinal mode there exists a critical energy threshold El under which the solution of (34) is 
stable while for larger energies the solution may be unstable: we also show that different initial data 
with the same total energy give the same stability response. Not only this enables us to numeric ally 
compute the threshold El and to evaluate the power of the sufficient condition given in Theorems 
and 4.5, but also to define a flutter energy for each longitudinal mode as a threshold of stability. 


4.4 


Definition 4.6. We call flutter energy of the k-th longitudinal mode Tp k (that is, the solution of (31 )) 
the positive number E k being the supremum of the energies El such that the trivial solution of (34) is 
stable for both l = 1 and 1 = 2. 


Concerning the bridge model, Theorems 4.4 and |4.5| lead to the conclusion that 
if the internal energy E is smaller than the flutter energy then small initial torsional oscillations re¬ 
main small for all time t > 0, whereas if E is larger than the flutter energy (that is, the longitudinal 
oscillations are initially large) then small torsional oscillations may suddenly become wider. 


5. Numerical computation of the flutter energy 

First, by using the eigenvalues found in Table [lj we numerically compute the parameters ak and bk 
in ( |30| ). It turns out that all the ak are equal to 0.1 up to an error of less than 10~ 3 : for this reason, 
in Table[3]we quote the values of 10 4 (afc — 0.1). Moreover, all the bk are equal to 1.1 up to an error of 
less than 10 -1 : we quote the values of 10 2 (&fc — 1.1). 


k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

10 4 (o fc -0.1) 

0.05 

0.2 

0.45 

0.8 

1.24 

1.78 

2.42 

3.14 

3.96 

4.86 

5.85 

6.92 

8.06 

9.28 

10 2 (6 fc -l.l) 

4 

4.03 

4.09 

4.17 

4.27 

4.39 

4.54 

4.7 

4.89 

5.1 

5.32 

5.57 

5.83 

6.11 


Table 3. Numerical values of ak and bk- 


Then we compute ai and d^k as defined in (36). We find that both a\ 
di t k are as reported in Table [4j 


0.27, 02 ~ 0.27, while the 


k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

di,k 

9.27 

6.18 

6.18 

6.18 

6.19 

6.19 

6.19 

6.2 

6.2 

6.21 

6.21 

6.22 

6.23 

6.24 

d2,k 

6.18 

9.27 

6.18 

6.18 

6.19 

6.19 

6.19 

6.2 

6.2 

6.21 

6.21 

6.22 

6.23 

6.24 


Table 4. Numerical values of d\ p and d, 2 ,k- 
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In fact, 0 < k — d 2 ,fc ~ 10 4 (for k = 3,..., 14) so that, in Table 0 one does not see any differ¬ 
ence between these coefficients: we put however the “exact” numerical values in the below numerical 
experiments. 

Concerning the structural parameters and the value of 7 we first notice that the total length of the 
composed spring depicted in Figure [4] equals the height of the towers over the roadway, which is 72 m, 
see {2j. We denote by h the sag of the parabola drawn by the cables: since the shortest hangers at 
midspan were of 2 m we have h = 70 m. We denote by 77 and V respectively the horizontal and vertical 
components of the tension of the couple of cables; in the equilibrium position under the action of the 
weight of the cables, of the hangers and of the deck, the tension is given by the vector 


(H,V) = \(H,V)\!> = 


\(H,V)\ 


Ah 

l,-(L-2x) 


1 + \jj( l - 2 a 0 ] 

where T" denotes the tangent unit vector to the parabola. From this we deduce that 


W,v)\ 


\A + [t?( l - 2x )\ 


= H, 


V = 


\(H,V)\ 


\A + [t?( l - 2x )\ 


5 ^(L-2z) = ffg(L-2x). 


The horizontal component of the tension 77 can be supposed constant with respect to x and hence we 
may write 

4 h 

V(x) = 77-—r (L — 2x) for any x E (0, L). 

Therefore, the action of the cables per unit of length is given by 


V'(x) 



for any x E (0, L ). 


For more details on the behavior of the cables and how to derive the above formulas we address the 
interested readers to [23., Chapter 3]. 

If we consider a configuration corresponding to a displacement w from the equilibrium configuration, 
the increment of the action of the cables per unit of length is given by jjW- Taking into account that 
the weight of the bridge per unit of length is approximatively 83 kN/m , we have that 77 = 1.08 • 10 8 TV. 

We have seen in Section [ 2 ] that the coupled action of cables and hangers is nonlinear. Following ([!]) 
we suggest as possible force the quantity 


8 77 

U 


(w + w 3 ). 


Next we compute the action of the cables and of the hangers per unit of surface: since the hangers 
act on the region lx defined in <§ we obtain 


1500 8 H, 600077 

Taking into account that m = 635 kg/m 2 we come to the equation 


m wtt + TA 2 w + 


600077 

~U~ 


T (y)(w + w 3 ) = 0 


where T denotes the rigidity of the plate. 

In order to provide a reasonable value for the rigidity compatible with the parameters of the TNB, 
we start by considering the rigidity of the deck in the beam model. The rigidity of a beam is given 
by El where E is the Young modulus and 7 is the moment of inertia of the cross section with respect 
to the horizontal axis orthogonal to the axis of the beam and containing its barycenter; from [ 2 j we 
know that E = 2.1 • 10 11 Pa and 7 = 0.1528 m 4 . Then, the plate equivalent to the beam has a rigidity 
given by T = 2 t(i-o- 2 ) = 2 -937 • 10 9 Pa ■ m 3 , see [16] for more details on the comparison between 

the two models. From [16] we also recall that the rigidity of a plate is given by 12 (^'-o- 2 ) w ^ ere d 
is the thickness of the plate. In our case we can recover d from the value of F: indeed we have 
d = [12(1 — a^T/E] 1 / 3 = 0.544 m. We observe that this is not the thickness of the roadway slab, 
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which is known to be of 13 cm, see [2], p.13]. However, since the roadway slab was reinforced by a 
stiffening girder with an H-shaped section, the deck may be considered as a plate of thickness four 
times as much: d ~ 52 cm. Finally, we recall that a rectangular plate has to be considered as a 
rectangular parallelepiped of thickness d made of a homogeneous material with constant density and 
isotropic behavior. 

We have so determined the constants in the model problem ([9]) and, defining u as in (12), we come 
to the adimensional problem (14) where 


(40) 


vr 4 r 


7 = 


6000 HL 


= 5.17- 10~ 4 . 


Since we aim to show that the value of 7 plays a secondary role, we quote below the results for three 
different 7. 


Once all these parameters are fixed, we solve (32) for (p k ( 0) = A > 0 and y>l(0) = 0 for different 
values of A. E ach value of A yields the /c-th longitudinal mode 1 p k = Tp k at energy E(A, 0) > 0 
Definition 


see 


4.2 


In turn, we use <p k to compute the function Ai tk (t ) defined in (35) and we replace it 


into (34). We start with A = 0 and we increase it until the trivial solution £0 = 0 of (34) becomes 


unstable. Since this is a very delicate point, let us explain with great precision how we obtain the two 
sets of critical values of A (thresholds of instability) that we denote by A\(k) and A 2 (k). 


For a particular second order Hamiltonian system of the kind of (28), it is shown in [9] that there 
exist two increasing and divergent sequences {Ay }^ =0 (/ = 1,2) such that: 

• if A 6 S := U k{A% k , A 2k+1 ) then £0 is stable; 

• if A € U := U k (A 2k+1 , A 2k+2 ) then £0 is unstable. 

Moreover, the instability becomes more evident if A is far from S: in particular, if A E (A 2k+1 , A 2k+2 ) 
for some k ^ 0 and A 2k+2 — A 2k+1 > 0 is small, then it is hard to detect the instability of £o- 
Due to the already mentioned unpredictable behavior of the stability regions for general Hill equations 
(see mm), the same results appear difficult to reach for the particular Hamiltonian system ( |28[ ). It 
is however reasonable to expect that somehow similar results hold. In particular, we expect £0 to 
be “weakly unstable” whenever A belongs to a narrow interval of instability. By this we mean that 


nontrivial solutions of (34) blow up slowly in time. In other words, if A belongs to a narrow instability 
interval, then only a small amount of energy is transferred from the longitudinal mode Tpt to a torsional 
mode. From a physical point of view, this kind of instability is irrelevant, both because it has low 
probability to occur and because, even if it occurs, the torsional mode remains fairly small. In turn, 
from a mechanical point of view, we know from Scanlan-Tomko m that small torsional oscillations 


are harmless and the bridge would remain safe. For this reason, we compute the two sets of critical 
values Ai{k ) (l = 1,2) as the infimum of the first interval of instability having at least amplitude 0.2. 
These critical values measure the least height of the longitudinal mode Tp k which gives rise to a “strong” 
instability. 

Let us explain what we saw numerically. For 7 = 10 -4 we discuss here the stability of the lAth 
longitudinal mode 1^(4 for A 6 {0.79; 0.8; 0.81} with respect to the first torsional mode. In the pictures 
of Figure [7] we represent the corresponding solution of (34) with £(0) = £ 7 (0) = 1. 



Figure 7. Solutions of (34) with £(0) = £'(0) = 1, k = 14, l = 1, A G {0.79; 0.8; 0.81}. 
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When A < 0.75 the behavior of £ appears periodic with, basically, oscillations of constant amplitude. 
When A = 0.79 (left picture) the function £ still appears periodic but with oscillations of variable 
amplitude; this behavior always turned out to be the foreplay of instability. And, indeed, when A = 0.8 
we see (middle picture) that £(f) oscillates with increasing amplitude and reaches a magnitude of the 
order of 10 4 for t = 200. The same phenomenon is accentuated for A = 0.81 (right picture) where 
£(200) has an order of magnitude of 10 s . By increasing further A the magnitude was also increasing. 
This is how we found Ai(14) = 0.8. 

From the values of A\(k) and A 2 (k) one can compute the corresponding energies E\ and E 2 as given 
by (f33l): 


(41) 


E t (k) = (7MM + a k (l = 1, 2) 


Then we can compute the flutter energy of the fc-th longitudinal mode ip k (see Definition 4.6) as 


E k = mm{Ei(k),E 2 (k)} . 


In Table [5] we quote our numerical results for some values of 7. In the middle table we use (40), that 
is, the value of 7 obtained from the parameters of the TNB. We found completely similar behaviors 
also for 7 = 10~ 2 and 7 = 5- 10 -5 . 



k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 


A\(k) 

4.62 

2.96 

2.93 

2.85 

2.67 

2.31 

1.42 

>20 

>20 

>20 

0.89 

1.51 

2.02 

2.51 

A 2 (k) 

5.93 

9.23 

5.91 

5.87 

5.79 

5.63 

5.36 

4.92 

4.19 

2.62 

>20 

>20 

>20 

>20 

k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

Ai(fc) 

3.32 

2.12 

2.10 

2.04 

1.92 

1.65 

1.01 

> 20 

> 20 

> 20 

0.62 

1.08 

1.46 

1.82 

A 2 (k) 

4.26 

2.46 

4.25 

4.22 

4.15 

4.05 

3.86 

3.54 

3.01 

1.87 

> 20 

> 20 

> 20 

> 20 


k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 


A\(k) 

1.44 

0.91 

0.9 

0.88 

0.82 

0.69 

>10 

>10 

>10 

>10 

0.2 

0.44 

0.63 

0.8 

A 2 (k) 

1.87 

2.91 

1.86 

1.85 

1.82 

1.77 

1.69 

1.54 

1.3 

0.76 

>10 

>10 

>10 

>10 


Table 5. Numerical values of A\{k) and A 2 {k) when 7 = 10 4 (top), 7 = 5.17 • 10 
(middle), 7 = 10“ 4 (bottom). 


1-4 


Remark 5.1. From Table [5] we deduce that, if 7 = 10 3 or 7 = 5.17 • 10 4 then A\(k) > A- 2 {k) 
provided that k = 8,9, 10. If 7 = 10 -4 this happens provided that k = 7, 8,9,10. In these cases, 
the energy transfer occurs on the second torsional mode. On the contrary, for lower k, we have that 
A\(k) < A 2 (k). 


These results deserve several comments. First of all, we notice that the values of Ai(l) and A 2 {2) are 
somehow out of the pattern because they are strongly influenced by the definition of the coefficients 
d^ k in (36) which is fairly different if l = k and l 7^ k. If we drop the case k = l, from Table [5] we 
see that the maps k i-a Ai(k ) are strictly decreasing until some k where A[(k) becomes very large. In 
particular, when l = 2 the least amplitude of oscillation A 2 {k) (threshold of instability) is obtained for 
k = 10. This behavior is obtained for different values of 7. As pointed out in the Introduction, on 
the day of the TNB collapse the motions were considerably less than had occurred many times before. 
Table [5] gives an answer to question (Q3): the tenth longitudinal mode seems more prone to generate 
the second torsional mode. 
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6. Validation of the results 


6.1. Validation of the linearization procedure. Fix some k E {1,..., 14} and some l E {1,2}. 
The energies Ei(k) in (41) are computed in Section [5] with the following algorithm. Firstly we solve 
(32) and find the fc-th longitudinal mode Tp k at energy E{(f >Then we use Tp k to determine the 
function Ai k (t) in (35). Finally, we study the stability of the trivial solution of (34): if it is stable then 
E((f >< Ei(k) while if it is unstable, then E(<p > Ei(k). With different choices of the initial 


data (r/>Q, <j>\) we are so able to find both upper and lower bounds which can be as close as wanted. 
Overall, this algorithm means that we are solving the following system 


Vfc(*) + (7A*M + ak)<Pk{t ) + h<Pk(t) = 0 Vf > 0 

€"(t) + (7^,2 + a t + = 0 Vt> 0 

^fc(O) = </>o , <p' k (0) = (/>i 

.6(0) =£0, €' l (0)=£l 


where |ei| + |£21 "C |6ol + | dq | • The system (42) is obtained by putting to 0 all the pj with j / k and 
Tj with j ^ l in (28) and by linearizing the so found 2x2 system of ODE’s in a neighborhood of the 
solution (jp k . 0). 

One may wonder if this method to determine Ei(k) does not suffer from this linearization and give 
incorrect results. We have so tried a nonlinear approach and considered systems such as 


Vfc(6 + (7Mfc,i + o,k)p k {t) + bkvl(t) + ai£?(t)(pk{t) + a 2 £i(t)ipl(t) = 0 \/t > 0 

6"(6 + (7*4,2 + ai + + ft <pk(t)^i(t) + fo£ii(t) = 0 Vt > 0 

^fc(O) = , ¥4(0) = ^ 1 

.6(0) = £0) 6(0) = £ i; 


note that the equations in ( |43| ) differ from the ones in (42) by third order homogeneous polynomials 
with respect to and 6• We numerically solved (43) for different choices of the constants on and /% 
but, regardless of the choice, we always found that 


the solution £/ of (42) remains small for all t > 0 if and only if the solution £/ of (43) 

remains small for all t > 0. 


This behavior had to be expected in view of the theoretical result in |2Bj. 

This means that the computation of E[(k) (and of the flutter energy of the fc-th longitudinal mode) 
does not depend on the linearization procedure. Moreover, this shows that Definition |4.3| well charac¬ 
terizes the stability of the /c-th longitudinal mode with respect to the Z-th torsional mode. 


6.2. Different nonlinearities: a system of Hill equations. For any longitudinal mode k, the just 


described procedure with the specific nonlinearity chosen, see (20), yields the system (34) of uncoupled 


Hill equations (1 = 1,2). This enables us to study separately the stability of each of the two equations 


in (34) and to compute the thresholds Ai(k) as described above. 

let 


For different nonlinearities, other than (|20j) , it may happen that the corresponding equations in (34) 
(l = 1,2) remain coupled and give rise to a Hill system of the form 


(44) 


^"(t) + Ai t k(t)z*(t) — 0 


where S = (6,^2) an d Ai k is now a 2 x 2 periodic matrix. In this case, one should investigate the 
stability of the trivial solution E = (0,0) of (44). By [35l Theorem II p.270, vol.l] we know that 
this trivial solution is stable provided that the trivial solutions of a family of related Hill equations 
are stable. Therefore, even for different nonlinearities or for a larger number of torsional modes, the 
stability of the /c-th longitudinal mode with respect to the torsional components may be performed 
through a finite number of stability analysis of some Hill equations. 
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6.3. The coupled nonlinear system. In this section we discuss the nonlinear system (28) with 
7 = 5.17 • 10~ 4 as for the TNB, see (40). We provide an alternative notion of stability for a solution of 


(28) having the torsional components Ty and t 2 identically equal to zero. Then we discuss the stability 


of the longitudinal components. Our purpose is to explain to which extent the decoupling procedure 


applied in Section 4.2 provides an accurate description of the phenomena. 

For our convenience we slightly modify some notations. We denote by wy,... ,wy± the first fourteen 


eigenfunctions (see Section 4.2) which are even with respect to y, i.e. Wk = Wk for any index k G 
{1,..., 10} and w n = wy 2 , w 12 = 1 x 7 . 3 , wy 3 = wyy, w 14 = 1 x 715 . Then we denote by w 15 = u>io and 
wiq = w i6 as the first two eigenfunctions which are odd with respect to y. All the eigenfunctions are 
normalized in L 2 (fi). Similarly we relabel the functions ipy,, ip 14 and Ty,T 2 as flk for k G {1,..., 16} 
and the corresponding eigenvalues ft k for any k G {1,, 16}. Then, for T as in (20), we introduce the 
constants 

Ay. = I T (y)wl dxdy for any k G {1,..., 16} 

Jn 

and for any jy,j 2 ,j 3 16}, jy ^ 32 > J3 and k G (1,..., 16} 


In T (.y) w ji u ’j 2 w j 3 w k dxdy 


(45) 


B 


'hhhk 


= 


if 3 1 = 3 2 = j’.i 


3 In T (y) w ji u ’j 2 w j 3 w k dxdy if jy > j 2 = 33 or j\ = j 2 > j 3 


. 6 In T (y) w h w j 2 w j 3 w k dxdy if jy > j 2 > j 3 . 

With these notations we may write ( |28[ ) and the corresponding initial conditions in the form 

16 

+ (7 idk + Af.)<p k (t) + B jij2j 3 k = 0 , k G {1,..., 16} , 


(46) 


jl,32 ,j 3 = 1 


<p k (0) = (f>^, <p' k (0) = 4>i , k G {1,..., 16} . 

We observe that if k = 15,16 then Bj 1 j 2 j 3k = 0 for any ji, j 2 , j 3 G (1,..., 14} such that j\ ^ j 2 ^ j 3; 


therefore if we choose 4>q 5 = = 4> o 6 = 4>\ h = 0) then the solution of (46) satisfies £75 = ip±Q = 0. 

This means that no torsional oscillation may appear if such a kind of oscillation equals zero at t = 0. 
We summarize this fact in the following 

Proposition 6.1. The unique solution (tpi,..., tpi^) of (46) with </>q 5 = cj)\ 5 = q 6 = (fi = 0 satisfies 
77 . 5 (f) = <Pw{t) = 0 for any f Gl. 


The situation is completely different if we do not assume that q i>g 5 = 7l’ 5 = = ^i 6 = in this case 

the first fourteen equations are coupled with the last two since the B k j 2 j 3k may be different from zero 
if k G {15,16} and j 2 ,j 3 ^ 14. This suggests to make precise how an oscillation mode may influence 
the others. 


Definition 6.2. Let us consider system (46). 

(i) We say that the k-th mode influences the j-th mode (j f k) if B kk kj f 0. 


(ii) 

(hi) 


We say that the k\-th and k 2 -th modes (k± > k 2 ) influence the j-th mode (j fL {k\,k 2 }) if at 
least one between B kl k 1 k 2 j or B k 1 k 2 k 2 j different from zero. 

We say that the ky-th, k 2 -th, k 3 -th modes (k± > k 2 > k 3 ) influence the j-th mode (j 0 
{ki, k 2 , k 3 }) if Bk 1 k 2 k 3 j 0. 


If some modes influence the j-th mode it may happen that, in system (46), <pj ^ 0 even if q, fpf) = 
(0, 0). We state a result which explains how longitudinal modes influence each other. 


Proposition 6.3. The following statements hold true: 

(i) Let k,j G {1,..., 14} with j k. If j = 3 k then the k-th mode influences the j-th mode; 

(ii) Let ky,k 2 ,j G {1,..., 14} with k\ > k 2 and j $ {k\. k 2 }. If j G {2ky + k 2 , 2ky — k 2 , ki + 2k- 2 , \ky — 
2k 2 \} then the ky-th and k 2 -th modes influence the j-th mode. 
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(iii) Let k\, k 2 , k 3 ,j € {1, ..., 14} with k\ > k 2 > k 3 and j 0 {ki, k 2 , k 3 }. If j E {h + k 2 + k 3 , k\ + 
k 2 — k 3l k\ — k 2 + k 3 ,\ki — k 2 — A&l} then the k\-th, k 2 -th and k 3 -th modes influence the j-th 
mode. 

Proof. In order to treat (i)-(ii) together we assume that k± ^ k 2 where also equality is admissible. 
Then, under the assumptions of the proposition we infer that either B^kx k^j 7 ^ 0 or B^k 2 k 2 j 0 0 since 
either sin(/cix) sin(Aqx) sin(^ 2 x) sin(jx) dx / 0 or sin(fcix) sin(/c 2 x) sin(/c 2 x) sin(jx) dx 0 thanks 
to Werner formulas; moreover Vk{y) > 0 for any y E ( — j-§o ’ ifo) anc ^ ^ ^ { 1 , - -., 14}, with Vk as in 
Section 4.2 The case (iii) can be treated in a similar way. □ 


Remark 6.4. Proposition | 6 .3| strongly depends on the nonlinearity h introduced in (20). That the k -th 
mode influences the j-th mode when j = 3k (statement (i)) depends on the fact that the nonlinearity 
h(y, u ) contains a cubic term. If the cubic term is replaced with a different term then the influences 


may vary. A similar discussion is valid for (ii)-(iii) in Proposition 6.3 

we consider the case where 


To illustrate Proposition 


6.3 


= 1 , 


= 0 and d>n = 


$ = 0 for 


all k E {2,..., 16}. According to Proposition 6.3, the first longitudinal mode influences only the third 
longitudinal mode but this does not mean that only 0i and 03 are nontrivial in (46). As one can see from 
Figure [ 8 ] below, where a numerical solution of (46) is obtained, also the functions 0s, 07 , 0g, 0n, 013 
are not identically equal to zero while all the other longitudinal modes and the two torsional modes are 
identically equal to zero. 



Figure 8. The solution of (46) with 0q = 1,0} = 0, (0q, 0 1) = (0) 0) f° r any A; E {2,..., 16}. 


Let us briefly discuss these results. From Proposition 6.3 (i) with k = 1 and j = 3, we deduce that 


01 influences tp 3 ; by Proposition 6.3 (ii) with k\ = 3 and k 2 = 1 we have that 0i and ip 3 influence <j%. 
Then 0i and 05 influence 07. In the next stages 0i and (97 influence 0g, while 0i and 0g influence 
011 ; finally 0 i and p \ \ influence 0 i 3 - 

As a second example we consider 0j} = 1, 0} = 0 and (0§, 0i) = (0,0) for any k E {1,..., 16} with 
k 0 2. A similar phenomenon occurs as one can see from Figure [9] below. One can observe that only 
the functions ip 2 , 06, 0io, 0 i 4 are not identically equal to zero. 

A further consequence of our approach deserves to be emphasized. We have truncated an infinite 
dimensional system up the least 16 modes: although this truncation is legitimate for energy reasons 
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~£3 

-<Pa 


0.4 

0.3 

0.2 


- 0.1 

- 0.2 

- 0.3 

- 0.4 


-<05 

~<P6 


- ify 

- y?8 




- £15 

~<P 16 


Figure 9. The solution of (|46j) with = 1, (f>\ = 0, {4>q,4>i) = (0,0) for any k G 

|l.16}, M 2. 

(see 0), it may lead to some small glitches. Since we are considering a truncated system with fourteen 
longitudinal modes and two torsional modes, it happens that if k E {5,..., 14} then the mode 
influences no other modes, hence the solution of (46) with ((/)§, (0,0) and (£ 0 , (/)■[) = (0,0) for 


any j G {1,..., 16}, j / k, satisfies <pk ^ 0 and t pj = 0 for any j G (1,..., 16}, j k. This is another 

consequence of Proposition |6.3| (i) which can be summarized in the following 


Proposition 6.5. Let k G {5,..., 14}. Suppose that (</>q, (j>\) (0, 0) and that (</>q, 0}) = (0,0) for any 

j G {1,..., 16}, j / k. Then the corresponding solution of ([46} satisfies 


<Pj(t) = 0 for any f G R and j G (1,..., 16}, j ^ k . 

6.4. Comparison of the results for the coupled and uncoupled systems. In this section we 


show that the responses of the truncated system (46) are in line with the theoretical and numerical 
responses obtained from the Hill equation (34) in Sections 4.3 and [5] 

We first provide the definition of stability of solutions of (46) satisfying (</>q 5 , <f\ 5 ) = (</> q 6 , = (0, 0), 

namely solutions of 


(Po) 


16 


Tk(t) + (7 hk + A k )<p k (t) + B hj2j 3 k Til ( t)<p j3 (t) = 0 , k G {1,..., 16} 

k G {!,•••, 14} 


31,32 ,j 3 = 1 


^fc(o) = 4 >q , ^(o) = ^, 

^15(0) = ^is(O) = <P16 (0) = Tw (0) = 0 • 


In the sequel a solution of ( Pq ) will be denoted by <ho = (<Pi, ■ ■ ■, 7 ’ 14 ,0, 0), see Proposition |6.1| We also 
denote by = {fp \,..., a general solution of (46). 


Definition 6.6. Let 


(Po)- 


t>Q ,..., </>q 4 and </>},..., (f\ 4 be fixed and let 4>o be the corresponding solution of 
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(i) We say that a solution 4>o of (Pq) is stable with respect to the first torsional mode if for any 
e > 0 there exists 5 > 0 such that if max{|^Q 5 |, |</>} 5 |} < 8 and (0q 6 ,<^>{ 6 ) = (0,0), then the 
solution of (46) corresponding to (f> q, ..., </>J 6 and (f\,, f>\ 6 satisfies 

\<pi 5 (t)\<e, |^ 15 (t)| < £ foranyteR. 


(ii) We say that a solution 4>o of (Pq) is stable with respect to the second torsional mode if for any 
e > 0 there exists 6 > 0 such that if max{|</>Q 6 |, 10i 6 1} < 5 and , (f>\ 5 ) = ( 0 , 0 ), then the 
solution of (46) corresponding to c/)q, ..., <f}^ and (j)\,..., 4> 1 l b satisfies 

|^i 6 (i)| < £, \<Pie(t)\<e foranyt&R. 


Definition 6.6 should be compared with Definition |4.3[ So far, to find sufficient conditions on <f>o 
which guarantee its stability with respect to one of the torsional modes (according to Definition 6 . 6 ) 
seems out of reach. However, Theorems 4.4 and 4.5 suggest the following 

Conjecture 6.7. Let l G {1,2}. There exists Ai> 0 such that if 

max{|<^o|,..., \&\..., |^i 4 |} < Ai 

then the corresponding solution 4>o of (|H)|) is stable with respect to the l-th torsional mode. 


Our purpose is to support Conjecture |6.7| from a numerical point of view under suitable restrictions on 
the values of the initial conditions 4>q, ..., 0 q 4 , ..., </>} 4 in (Po )• For more numerical results concerning 

Conjecture 6.7 we refer to [7j. 

For any k G {1,..., 14} we consider the problem 

(Po,k) 

16 

Vj(t) + (7 H + + ^2 Bj lj2j3 j fih(t)fi h {t)fi j3 (t) = 0 , j G {1,..., 16} , 


h ,32,33=1 
31>j2>33 


V? fc (0) = H, <Pfc(0) = 0 , 
Tj ( 0 ) = <p'j ( 0 ) = 0 , 


j g {l,..., 16} ,j k. 


Then for any k G {1 

r 


(P6,k,l 


, 14} and l G {15,16} we also consider the following perturbation of (Po,fc) 
16 

Tj(t) + (7 H B hhhj Th(t)Tj 2 {t)<Pj 3 {t) = 0 , j G {1,..., 16} , 


h ,32,33=1 
31>J2 


ip k (0) = A 1 (p' k ( 0) = 0, 
£t(°) = $(°) = 

l^-(O) = ^'( 0 ) = 0, 


j e {!,••■, 16 },j fik,j fil. 


In the particular case of solutions of (To.fc), Conjecture 6.7 reads 

for any l G {1, 2} and k G {1,..., 14} there exists Afik) > 0 such that for any A G (0, Afik)) the 
solution <3?o of ( Pqj- ) is stable with respect to the l-th torsional mode . 


Numerical simulations give the approximate values of Afik) as reported in Table | 6 j 

Remark 6.8. From Table [ g] we see that A\{k) > A 2 (k) provided that k = 8,9,10 whereas for lower k 
we have that A\{k) < ^(L). 

Let us explain how we obtained Table [ 6 j As a particular example we take the value of Ai(5). In 
Figure 10 we plot the graph of the component ^5 of the solution 4 >q of (Po,k I with k = 5 and A = 1.92. 
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k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

A\(k) 

2.38 

1.89 

3.83 

2.27 

1.92 

1.66 

1.02 

> 10 

> 10 

> 10 

0.62 

1.08 

1.46 

1.83 

Mk) 

5.17 

4.38 

4.94 

8.08 

4.18 

4.05 

3.87 

3.59 

3.10 

1.87 

> 10 

> 10 

> 10 

> 10 


Table 6 . Numerical values of A\(k) and A 2 {k) when 7 = 5.17 • 10 4 . 



Figure 10 . The component <^5 of the solution Tq of (Po,k) with k = 5 and A = 1.92 


The instability phenomenon appearing for A = 1.92 is represented in Figure [Tl] where we see that 
the amplitude of the oscillations of £75 is fairly wide for large t even if we choose a small value of 5 in 


(Ps,k,i)- 









Figure 11 . The components of the solution <f> of (Ps t k,i) with k = 5, l = 15, A = 1.92 
and 6 = 5 • 10 -4 . 


On the other hand, if we choose A = 1.91 we see in Figure [l2| that the amplitude of the oscillations of 
£>15 remains of the same order of magnitude also when t is large. This is how we obtained Ai(5) = 1.92 
in Table |U 

The comparison between Table [b] and Table [ 5 ] (middle table for 7 = 5.17 • 10 -4 ) deserves several 
comments. The values of A\(k) when k G {5,.... 14} \ {8, 9,10} and of A-^ik) when k G {5,..., 10} are 
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Figure 12 . The function <^15 corresponding to the solution of (Ps,k,i) with k = 5, 
1 = 15, A = 1.91 and <5 = 5- 10 " 4 . 


very close to the corresponding ones of A\{k) and A 2 {k) in Table [5j The reason is that for these values 
of k the solution 4>o of (Po,fc) satisfies = 0 for any j 7 ^ k , see Proposition 6.5 hence the equations 
are decoupled. 

Then, with the same choice of k, fix S > 0, l £ {15,16} and let be the corresponding solution of 
( Ps,k,i )• If 5 is small enough the function is expected to be close to 4>o until <pi remains small together 
with its derivative. And indeed, if we chose k = 5, l = 15, A = 1.92, <5 = 5-10 4 , we numerically obtain 


the plots of Figure 11 


We also observe that for larger values of A, say A = 1.94, when the fifteenth component of is 
no longer “negligible” the values of the fifth component p 5 of <£,5 are considerably different from the 
values of the fifth component of 4>o, at least for large t: in the latter case ^5 is periodic while in the 
former case £>5 exhibits a variation in the amplitude of the oscillations soon after the amplitude of the 


oscillations of the fifteenth component has become relatively large, as one can see from Figure 13 



Figure 13. The functions tp 5 and ip 15 corresponding to the solution <f> of (Po,fc) with 
k = 5 and A = 1.94. 


The critical values Ai( 8 ), Ai(9), Ai(10), ^(ll), ^ 2 ( 12 ), ^2(13), ^2(14) are very large in both Tables 
[5] and [ 6 j No evident instability behavior was detected when considering solutions of ([32} and ([34} on 


one hand and solutions of (Ps,k,l) on the other hand. 

It remains to discuss the cases where k 6 {1, 2, 3,4} for which the values of A\(k) and A 2 {k) in Table 
[ 6 ] are quite different from the values of A\(k) and A 2 {k) in the middle Table 5j By looking at Figures [ 8 ] 
and [9] one can see the graphs of the components of 4>o, the solution of (Po,fc)> respectively in the cases 
k = 1 and k = 2 (for the two remaining cases k = 3,4 we refer to [7J). In the two figures one observes 
that the k -th component of 4>o is not the only large one, there are several other large components. 
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In this case, the reduction of (46) to equations (32) and (34) is not a good approximation of the real 
situation. 


6.5. Thresholds of stability in the TNB. In this subsection, starting from the values obtained in 
Table [6j we provide estimates on the thresholds of stability for the TNB by using the isolated system 
(14) with the values of the parameters given in Section[5j According to the notations of Subsection |6.3[ 
we denote by w k , k G {0,..., 14} the first fourteen longitudinal eigenfunctions normalized in L 2 (fl). 
For any k G {0,..., 14} and l G {1,2} we choose the initial conditions in ( |14| ) in the form 

(47) u 0 (x,y) = Ai(k)w k {x,y ), u 1 (x,y) = 0 for any (x,y) G (0,7r) x (—£,£) 


where the Ai(k) are as in Table [bj Our purpose is to provide, for any k G {1,..., 14}, the values 
measured in meters of the stability threshold for the energy transfer from the fc-th longitudinal mode 
to a torsional mode. This will give an idea of the initial displacement of the deck sufficient to activate 
the torsional oscillations. To this end, we first give in Table [7] the approximate values of the L°°-norms 
of the L 2 -normalized eigenfunctions. 


k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

\\ixk ||l°° 

3.897 

3.899 

3.899 

3.900 

3.901 

3.902 

3.904 

3.905 

3.907 

3.909 

3.912 

3.914 

3.917 

3.920 


Table 7. Approximate values of ||ii>fc||L°° for k G {1,..., 14}. 


By (47), Table [gJ Table [t] and (12) with k\ = k ‘2 = eo( ^ H we obtain the following table which shows 
that the instability amplitude has the same order of magnitude as observed for the TNB, see 


k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

IKHl- 

9.27 

7.37 

14.93 

8.85 

7.45 

6.48 

3.98 

> 40 

> 40 

> 40 

2.43 

4.23 

5.72 

7.17 

IKHl- 

20.15 

17.08 

19.26 

31.51 

16.30 

15.80 

15.11 

14.02 

12.11 

7.31 

> 40 

> 40 

> 40 

> 40 


Table 8. The L°°-norm of uq measured in meters corresponding to the stability thresh¬ 
old of the fc-th longitudinal mode with respect to the first torsional mode (first line) and 
the second torsional mode (second line). 


7. Proof of Theorem 14.41 

Assume that 1 ^ k ^ 14 and l G {1,2} are given. We introduce two constants which will be useful in 
the sequel. We set 

(48) Si := 7^2 + d t , p k ■= 7Mfc,i + • 

Using the definition of pk in (48), for any E > 0 we put 

sj p 2 k + 4 b k E ± p k 


(49) 


A k ±(E) = 


bk 


Then ([33]) reads 
(50) 


(A) 2 = j(K( E ) + vlW-(E)~vl) 


Hence, since any k -th longitudinal mode <p k satisfies (33), we deduce 
(51) 


= ^A 1(E). 


Since in the equation (32) there are only odd terms, the period Tk(E ) of tp k is the double of the 
width of an interval of monotonicity of tp k . As the problem is autonomous, we may assume that 
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and ^.(0) = 0. By symmetry and periodicity we then have that (p k (T k / 2) = 


V>k( 0) = - 

and Tp' k {T k / 2) = 0. In the first interval of monotonicity of Tp k we may take the square root of (50) and 
obtain 


(52) 


¥>'*(*) = \/y ( A +(Y> + <4W)( A -Y) - <Pfc(*)) 


Vi £ 


T k 

°’Y 


By separating variables, and integrating over the interval (0,T k /2) we obtain 


T k (E) 


Hl^fcllc 


ds 


h J -iis 


llvjc 


^(A*(.E) + s 2 )(A*(£)-s 2 ) 

Using the fact that the integrand is even with respect to s and with a change of variable we get 

4\/2 f 1 ds 


(53) 


In particular, 


T k (E) = 


■'0 ^(A* (£) + A* (£)s 2 )(l - S 2 ) 


(54) 


27T 

the map E i—)■ T k (E) is strictly decreasing and lim T k (E ) = T k ( 0) =- 


Let us prove the following sufficient condition for the torsional stability of ip k . 
Lemma 7.1. If there exists an integer rn such that 


(55) 


4m 2 7r 2 


4 (m + l) 2 7r 2 


< T k (E) 2 < 


' ^d/ + dz, fc A-Y) 

f/ien is stable with respect to the l-th torsional mode. 


Proof. With the initial conditions £(0) = £(0) = 0, the unique solution of (34) is £ = 0. The statement 
follows if we prove that (55) is a sufficient condition for the trivial solution £ = 0 to be stable in the 
Lyapunov sense, namely if the solutions of (34) with small initial data |£(0)| and |£(0)| remain small 
for all t ^ 0. 

Since Tp k is TVperiodic, the function Tp\ is T k / 2-periodic. Then Ai )k {t) is a positive T^/2-periodic 
function and a stability criterion for the Hill equation due to Zhukovskii [3B], see also ESI Chapter 
VIII], states that the trivial solution of (34) is stable provided that 

Vt £ M 


4m 2 7r 2 . . . 4(m + l) 2 7r 2 


T k (Ef 


T k {Ef 


for some integer m ^ 0. By recalling (48) and the definition of A\ k in (35), this condition is equivalent 
to 

4m 2 7r 2 „ 4(m + 1) 2 7 t 2 


T k (Ef 


^ hi + d kk (p k (t) ^ 


T k (Ef 


Vt £ 


In turn, by invoking (51) we see that the latter is equivalent to (55). 


□ 


Remark 7.2. For E = 0 the right hand side of (55) equals 77 which is strictly larger than the 

left hand side of (55). But the right hand side of (55) is strictly decreasing with respect to E and tends 
to 0 as E — y oo. Therefore, the interval defined by (55) is empty for sufficiently large E. 


Let 

(56) 


m = max < /c £ N; k < \ — 

Pk 


so that, by (54), 4n V' 71 " J < T k { 0) 2 . We then infer that there exists E\(l, k) > 0 such that 


4m 2 7r 2 


< T k (E ) 2 VE^Ei&k). 


(57) 









































27 


This gives a sufficient condition for the first inequality in (55) to be satisfied. 
In view of (56) and of the assumption (|37j) we know that 


4-7T‘ 


T k (0) 2 = - < 


4(m + l) 2 7r 2 4(m + 1) 2 7 t 2 


Pk A* (0) 

Then the continuity of the maps E e -> T k (E) and E > A l f(E) implies that there exists E 2 (l,k) > 0 
such that 

4(m + 1) 2 7 t 2 


(58) 


T k (Ey ^ 


Si + d^ k A k _{E) 


VE^E 2 (l,k) 


This gives a sufficient condition for the second inequality in (55) to be satisfied. 

We may now conclude the proof of Theorem 4.4. We choose m as in (56), we define E\{l,k ) and 
E 2 (l,k) as in (57)-(58), and we take 

E^E l k := min{£i(Z,fc),£ 2 (Z,fc)} 


so that both (57) and (58) are satisfied. Then Lemma 7.1 tells us that Tp k is stable with respect to the 
Z-th torsional mode. By taking A = A k _ as in (49) we readily obtain the L°°-bound for tp k and the proof 
of Theorem |4.4| is so completed. 


The just completed proof also enables us to give the following quantitative version of Theorem 4.4 


Theorem 7.3. Fix 1 ^ k ^ 14, l € {1, 2} and let E l k := min{Td(Z, k),E 2 (l, k)}, see (57)-(58). Then the 
k-th longitudinal mode Tp k at energy E((f>is stable with respect to the l-th torsional mode provided 
that E ^ E l k or, equivalently, provided that ^ A k _(E l k ) where A k i is defined in (49). 

8. Proof of Theorem EH 

The proof of Theorem 4.5 follows the same lines of the proof of Theorem |4.4[ see Section [7J But the 
continuity of the maps involved is not enough to obtain the desired inequality and a different stability 
criterion is needed. 

Note that by (48) the condition (38) reads 
(59) 


Si 


— = (m + l ) 2 , 
Pk 


The right hand side of (|53|) is an elliptic integral of the first kind: after the further change of variables 

d<j) 4 


s = cos it may be written as 

4^2 r ' 2 


(60) T k (E) = 


where 


t/2 


dcj) 


'/h Jo ^A k + (E) + A k l(E) cos 2 yj p\ + 4:b k E •'° V 1 ~ Pk(E) sin 2 </> 


(61) 


Pk 


Pk(E) = - 1 - - - 

y/A + MkE, 


°'2 


This enables us to compute the derivative of T k (E) for E = 0. 

Lemma 8.1. Let T k = T k (E) be the function in ( |60[ ). Then 

3 7T b k 

T t(° =- 5 /| ' 

2 Pk 


Proof. An asymptotic expansion of p k in (61) shows that 

Pk(E) ~ E as E 

Pk 
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Then we may also expand T k (E) in (60) as follows 


T k (E) 


i -He 


r 71-/2 


dcf> 


V~Pk \ Pk J J o , /l — %-E sin 2 

Pk 


4 ( 

2 vr / 


v 1 p\ E 

h 


(*7r/2 { L 

1 H- \e sin 2 4> ) d(f> 

^Pk 


\[Pk V Pk 


b k 


1 — —q .£/ ) I 1 H—-—o E 


4p 


This proves that 


and the statement follows. 


T t (E)=T t (0)- : PP^E + c(E) 

2 pf 


as E —> 0 . 


as E —> 0 


□ 


Due to the presence of odd terms, the solution of (32) has several symmetry properties that we 
summarize in the next statement. We also provide a pointwise upper bound for Tp k . 

Lemma 8.2. Let if be the unique (periodic) solution of the problem 

jV'W + p k if(t) + b k if 3 (t) = 0 Vt > 0 


(62) 


v>(0) = 0 , V>'(0) = V2E 


and denote by T its period. Then: 

(i) i />(T/4) = max* if{t), if(T/A — t) = if(T /4 + t ) for all t, if(T /2 + t) = —if(t) for all t; 

(ii) the following estimate holds 


f 2E /-1 T 

0 ^ if(t) ^ min J — sin (y/pkt ), \JAt(E) > V 0 ^ t ^ — . 


Proof. The symmetry properties in (i) are well-known calculus properties. 

In order to prove (ii), we observe that by conservation of energy and arguing as for (52), we find 

if 1 = ^2 E - p k if 2 - y t /; 4 VO ^t^j. 

By separating variables we then obtain 

d u 
7 = at- 

yJ2E - p k if 2 - 

in turn, by integrating over [ 0 , t] and dropping the fourth power term, we get 


1 


p/Pk 

This yields 


■ (rp* 1 r m di ^ 

arcsm t —m = . / — , = V / —t 

VV 2S ) pj2E Jo Jo 


di 


< 2 E- Pk e-k^ 


l 


= dr = t. 


2 E T 

i>{t) ^ \ —sin (y/~Pkt) - 

VPk 4 


while from (51) we know that if(t) ^ yjAf_{E) for all t. By taking the minimum, we infer the desired 
upper bound in (ii). 


□ 


The previous lemmas enable us to prove the following technical result. 
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Lemma 8.3. Assume that (39) and (59) hold. There exists E 2 (l,k ) > 0 such that if E ^ E 2 (l,k) then 

r T k{E)/2 — l /maxf Ai k (t)\ 

Jo /VW‘ i(+ 2 l 0 g (mm l4 , i ( t ))« (m+1) ^ 

T (E) 

Proof. Recall that A^ k is a ^ -periodic function. Up to a time translation, we may assume that Tp k 


solves (62); then we estimate 
rT k (E)/2 ,- 

J y Ai^{t) dt = 

as E —y 0 ~ 


rT k (E)/ 4 

2 yjA hk {t)dt = 2 


, - fT k (E)/4 r - 

y dt = 2 J y s i + di ik Jpl(t) dt 

T k (E) +[ Tk(E)/ \m 2 {y/Vkt)dt 

2 PkV°i J o 

3 A , 2d l,kE fT t (E) sm(^nT k (E)/2)\ 

T" ( r ‘ (0) - + 7^5 V~8-VS-J 





T fe (£)/4 


by Lenuna 8.2 
(63) by Lemma 8.1 




By Lemma 8.1 we also infer that 

/ y^T fc (u) \ 


Sill 






= sin 


by (54) = sin I 7 r — 


Pk V$i 

. fy/f*T k (0) 37T b k ^ t f ^\ 
m \—2 - Tpt E + 0{E) ) 

^E + o(E))=o( 1) as .E —>• 0 . 

4 Pk J 


By plugging this information into (63) we end up with 

rTk{E)/2 A7 Tv ^ 37Tbk , d l,k T ki.O) c , 

« T( T ‘ (0) ~^ £ j + VA + " ( 1 


/o 


(64) 

O: 

sharp bounds 


i iHnK ( , *' Z' rfi,fc 3(m+l)6 fe 

by@ ~ ( ro + 1 )’ r + ^|(sVT- 2 -, 

On the other hand, by (51), we know that the periodic coefficient A) k in (34) satisfies the following 


E as E —> 0 . 


Si ^ A itk (t) ^ 6 t + di tk A k _(E) Vt^O 
where 5i is dehned in (48). In fact, A^ k has only two critical points in [0, T k (E)/2) and 


nun 

te[ 0 ,T fe (E)/ 2 ) 


Ai' k (t) = 5i , max A^ k {t) = 8i + d i)k A k _(E) . 
t£[0,T k (E)/2) 


In particular this means that 

■o / m ax, V(V =log / 1 + ja A . )' 

V nnn t Ai ik (t) J V °l / 

and, by (59) and by taking advantage of the explicit expression in (49), we obtain 

, /maxf Ai k (t)\ di k k d [)k _ 2 di k 

(65) log —- , ; ' ~ -z-A_(E) ~ 2 - 4 - E = - fr^E as E 0 . 

\min t A itk (t)) Si Sip k (m + l) 2 p\ 

By combining (64) with (65) we finally infer that 

fT fc (E)/ 2 


[ T k(E)/2 j- — l /max; Ai k (t ) 

/ y Aj fe (t) dt + - log — - ’ ;; 

Jo V ’ 2 B V min t Ai,*(t) 


The 


w n . / 7 rdj,fc rfi.fc 3vr(m + l) 6 fc \ E ^ 

^ (tyi + 1)7T + ( —- — + -- 7 T ta ■ J —9 + o(E) as E >• 0 . 

\2(m + 1) (m+1 ) 2 4 J p\ 


statement then follows from assumption (39). 


□ 
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Next, we remark that the function 

fT k (E )/2 

E\-+ / 

Jo 


,/FU it- hog ( '■ maxt 

V lMJ 2 B Vmint A, )fc (t) 


(with a minus sign before the logarithm!) is continuous and, for E = 0, it is equal to 2 ° \JSi = (m+l)ir 


in view of (54) and (59). Therefore, there exists E\(l,k) > 0 such that 
r T k( E )/2 r- —— 1 fmax t Ai k (t) 

Jo V ,W 2 * \min t A lfk (t) 

By putting E l k := min{i?i(£, k),E 2 (l, k)} and by combining this result with Lemma 8.3 we infer that 


^ mn \/E^Ei(l,k). 


rT k (E )/2 


mn A 


JAA) dt-- log ( m "' t'f? 

v ’ 2 *\mm t A ltk (t) 


r T k (E )/2 


SC 


sJa^AJ) dt +\ lo g f 


/max t A, k (t) 


^ (m + l)7r V.E < E\ 


k * 


/o ’ * \min t Ai jk (t ) 

Then a stability result by Burdina ESI (see also Test 3, p.703]) allows us to conclude that the 
trivial solution of (34) is stable. By Definition |4.3| this means that the k -th longitudinal mode Tp k at 
energy E ^ E l k is stable with respect to the l -th torsional mode. This completes the proof of Theorem 

run 

Remark 8.4. Instead of the Burdina stability criterion, one may try to use different criteria which 


may need assumptions other than (39). For instance, the Zhukovskii criterion (already used in Lemma 
7.1) needs the assumption that 2 di ik ^ (m + l) 2 bk] since it is more restrictive than (39), it seems that 


the Burdina criterion performs better in this situation. But there exist many other criteria, see 


and some of them could allow to relax further (39). And, perhaps, some criterion would allow to drop 


any assumption of this kind after estimating directly di >k and b k . 

9. Conclusions: our answers to questions (Q1)-(Q2)-(Q3) 

In Section 4.3 we have obtained stability results for a finite dimensional approximation of ( |14[ ). We 
analytically proved that if a longitudinal mode is oscillating with sufficiently small amplitude then it is 
stable, that is, it does not transfer energy to torsional modes. The numerical results in Section [5] show 
that if a longitudinal mode is oscillating with sufficiently large amplitude then it is unstable, that is, 
it transfers energy to a torsional mode. These results are numerically validated in Section |6.3| thanks 


to a more precise approximation of (14), see (46); the small discrepancies between the two approaches 


are justified in detail. For more numerical simulations corresponding to other values of k and l we refer 


to [?]. Overall, recalling Definition 4.3, these results enable us to give the following answer to question 

(Ql). 

Longitudinal oscillations suddenly transform into torsional oscillations because when the flutter energy 

is reached the longitudinal mode becomes unstable with respect to a torsional mode. 

A few days prior to the TNB collapse, the project engineer L.R. Durkee wrote a letter (see El P-28]) 
describing the oscillations which were so far observed at the TNB. He wrote: Altogether, seven different 
motions have been definitely identified on the main span of the bridge, and likewise duplicated on the 
model. These different wave actions consist of motions from the simplest, that of no nodes, to the most 
complex, that of seven modes. On the other hand, we have repeatedly recalled that the day of the 
collapse the motions, which a moment before had involved a number of waves (nine or ten) had shifted 
almost instantly to two. 

In Section[3]we found the explicit form of both the longitudinal and torsional modes. We also obtained 
accurate approximations of the corresponding eigenvalues when the TNB parameters are considered. 
We have analyzed the first and second torsional modes although we have explained in Figure [6] why the 
cables inhibit the appearance of the first mode. This is confirmed by recent studies in m- In Remarks 
5.1 and 6.8 we emphasized that A\(k) > ^(/c) provided that k = 8,9,10 while for lower k we have 
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that A\(k) < A 2 {k). According to the above reported letter of Durkee, the TNB never oscillated with 
k = 8,9,10 before the day of the collapse. 

The results in [9] tell us that the energy transfer occurs when the ratio between the torsional and 
longitudinal frequencies is small (close to 1). And in Remarks 5.1 and 6.8 we saw that the energy of 
the k-th longitudinal mode for k = 8,9,10 transfers earlier to the second torsional mode rather than to 
the first. 

Overall, these results allow us to give the following answer to question (Q2): 


Torsional oscillations appear with a node at midspan because the cables inhibit the appearance of 
the first mode and the longitudinal modes prior to the switch to torsional modes have an instability 
threshold with respect to the second torsional mode smaller than with respect to the first torsional 
mode. 


Invoking again the results in one reaches the conclusion that the critical amplitude of oscillation 
of a longitudinal mode depends on the ratio between the torsional and longitudinal frequencies. Table 
[I] shows that this ratio reaches its minimum for the 10th longitudinal mode. And precisely the 10th 
mode was observed the day of the collapse prior to the appearance of torsional oscillations. Therefore 
our answer to question (Q3) is as follows. 

There are longitudinal oscillations which are more prone to generate torsional oscillations; in the 
case of the TNB the most prone was the 10th longitudinal mode. 
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